7-6 Shortness of breath

Patient reported shortness of breath to day 28.

Authors
Affiliations
Rob Mahar

University of Melbourne

James Totterdell

University of Sydney

Published

July 12, 2022

Preamble

Load packages
library(tidyverse)
library(labelled)
library(kableExtra)
library(cmdstanr)
library(posterior)
library(bayestestR)
library(bayesplot)
library(matrixStats)
library(patchwork)
library(lubridate)
library(ggdist)

theme_set(theme_classic(base_size = 10, base_family = "Palatino") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank()))
bayesplot_theme_set(theme_set(theme_classic(base_size = 10, base_family = "Palatino") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank())))

color_scheme_set("red")
options(digits = 4)
Prepare dataset
devtools::load_all()
all_dat <- read_all_no_daily()

# Anticoagulation set (intention-to-treat, including missing)
acs_itt_dat <- all_dat %>% 
  filter_acs_itt() %>%
  transmute_model_cols_grp_aus_nz()

# Anticoagulation set (intention-to-treat, excluding missing)
acs_itt_nona_dat <- acs_itt_dat %>% 
  filter(!is.na(out_sob))  

# Concurrent enrolments for C2
acs_itt_concurc2_dat <- acs_itt_dat %>%
  filter_concurrent_intermediate()
acs_itt_concurc2_nona_dat <- acs_itt_concurc2_dat %>% 
  filter(!is.na(out_sob))

# Concurrent enrolments for C3
acs_itt_concurc3_dat <- acs_itt_dat %>%
  filter_concurrent_std_aspirin()
acs_itt_concurc3_nona_dat <- acs_itt_concurc3_dat %>% 
  filter(!is.na(out_sob))

# Concurrent enrolments for C4
acs_itt_concurc4_dat <- acs_itt_dat %>%
  filter_concurrent_therapeutic()
acs_itt_concurc4_nona_dat <- acs_itt_concurc4_dat %>% 
  filter(!is.na(out_sob))
Load the Stan models
model_full <- cmdstan_model(file.path("stan", "binary", "logistic_site_epoch.stan"))
model_site_only <- cmdstan_model(file.path("stan", "binary", "logistic_site.stan"))
model_epoch_only <- cmdstan_model(file.path("stan", "binary", "logistic_epoch.stan")) 
model_fixed_only <- cmdstan_model(file.path("stan", "binary", "logistic.stan"))
Helper functions
make_summary_table <- function(dat, format = "html") {
  tab <- dat %>%
    mutate(CAssignment = factor(
      CAssignment, 
      levels = c("C1", "C2", "C3", "C4"),
      labels = str_replace(intervention_labels()$CAssignment[-1], "<br>", " "))) %>%
    group_by(CAssignment) %>%
    summarise(
      Patients = n(),
      Known = sprintf(
        "%i (%.1f)", sum(!is.na(out_sob)), 100 * mean(!is.na(out_sob))),
      Missing = sprintf(
        "%i (%.1f)", sum(is.na(out_sob)), 100 * mean(is.na(out_sob))),
      `Shortness of breath day 28` = sprintf(
        "%i (%.1f)", sum(out_sob, na.rm = TRUE), 100 * mean(out_sob, na.rm = TRUE)),
    ) %>%
    bind_rows(
      dat  %>%
        group_by(CAssignment = "Overall") %>%
        summarise(
          Patients = n(),
          Known = sprintf(
            "%i (%.1f)", sum(!is.na(out_sob)), 100 * mean(!is.na(out_sob))),
          Missing = sprintf(
            "%i (%.1f)", sum(is.na(out_sob)), 100 * mean(is.na(out_sob))),
          `Shortness of breath day 28` = sprintf(
            "%i (%.1f)", sum(out_sob, na.rm = TRUE), 100 * mean(out_sob, na.rm = TRUE)),
      )
    ) %>%
    mutate(CAssignment = fct_inorder(CAssignment)) %>%
    gather(key, value, -CAssignment, factor_key = T) %>%
    spread(key, value)
  colnames(tab)[1] <- "n (\\%)"
  if(format == "latex") {
    colnames(tab) <- linebreak(colnames(tab), align = "c", linebreaker = "<br>")
  }
    kable(tab,
      format = format,
      align = "lrrrrr",
      escape = FALSE,
      booktabs = TRUE
    ) %>%
    kable_styling(font_size = 10, latex_options = "HOLD_position")  %>%
    row_spec(nrow(tab), bold = TRUE)
}

make_stan_data <- function(
    dat, 
    vars = NULL, 
    outcome = NULL, 
    beta_sd = NULL, 
    ctr = contr.orthonorm) {

  X <- make_X_design(dat, vars, ctr)
  nXassign <- sum(grepl("rand", colnames(X))) - 1
  beta_sd <- c(2.5, rep(1, nXassign), beta_sd)

  y <- dat[[outcome]]
  N <- dim(X)[1]
  K <- dim(X)[2]  
  
    epoch  <- dat$epoch
  M_epoch  <- max(dat$epoch)
    region <- dat[["ctry_num"]]
  M_region <- max(region)
    site <- dat[["site_num"]]
  M_site <- max(site)
  region_by_site <- region_by_site <- dat %>% 
    dplyr::count(ctry_num, site_num) %>% 
    pull(ctry_num)
  
  list(N = N, K = K, X = X, y = y,
       M_region = M_region, region = region,
       M_site = M_site, site = site,
       M_epoch = M_epoch, epoch = epoch,
       region_by_site = region_by_site,
       beta_sd = beta_sd)
}

make_stan_data_concurrent <- function(
    dat, 
    vars = NULL, 
    outcome = NULL, 
    beta_sd = NULL, 
    ctr = contr.orthonorm) {
  
  X <- make_X_design(dat, vars, ctr, includeA = FALSE)
  nXassign <- sum(grepl("rand", colnames(X))) - 1
  beta_sd <- c(2.5, rep(1, nXassign), beta_sd)
  y <- dat[[outcome]]
  N <- dim(X)[1]
  K <- dim(X)[2]  
  list(N = N, K = K, X = X, y = y, beta_sd = beta_sd)
}

Descriptive

Anticoagulation

Anticoagulation
save_tex_table(
  make_summary_table(acs_itt_dat, "latex"),
  "outcomes/secondary/7-6-anticoagulation-summary")
make_summary_table(acs_itt_dat)
n (\%) Patients Known Missing Shortness of breath day 28
Low dose 610 577 (94.6) 33 (5.4) 115 (19.9)
Intermediate dose 613 584 (95.3) 29 (4.7) 110 (18.8)
Low dose with aspirin 283 271 (95.8) 12 (4.2) 59 (21.8)
Therapeutic dose 50 44 (88.0) 6 (12.0) 11 (25.0)
Overall 1556 1476 (94.9) 80 (5.1) 295 (20.0)
Table 1: Summary of outcome by anticoagulation intervention.

Age

Relationship age to shortness of breath
agedat <- acs_itt_dat %>%
  dplyr::count(out_sob, AgeAtEntry) %>% 
  spread(out_sob, n, fill = 0) %>% 
  mutate(
    n = `0` + `1` + `<NA>`,
    p = `1` / (`1` + `0`))
agemod <- glm(
  cbind(`1`, `0`) ~ AgeAtEntry, 
  data = agedat, 
  family = binomial())
agedat <- agedat %>%
  mutate(
    ypred = predict(agemod, newdata = agedat, type = "response")
  )
p1 <- ggplot(agedat, aes(AgeAtEntry, n)) +
  geom_col(colour = "grey40", fill = "grey40") +
  geom_vline(xintercept = 60, linetype = 2) +
  labs(y = "Number of\nparticipants",
       x = "Age at entry")
p2 <- ggplot(agedat, aes(AgeAtEntry, p)) +
    geom_point() +
    geom_vline(xintercept = 60, linetype = 2) +
    geom_line(aes(y = ypred)) +
    labs(y = "Proportion\nshortness of breath\nday 28", x = "Age at entry")
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-6-age.pdf"), p, height = 2, width = 6)
p

Country

Relationship country to outcome
tdat <- acs_itt_dat %>%
  dplyr::count(Country = factor(
    country, 
    levels = c("IN", "AU", "NP", "NZ"),
    labels = c("India", "Australia", "Nepal", "New\nZealand")), out_sob) %>%
  group_by(Country) %>%
  spread(out_sob, n, fill = 0) %>%
  mutate(
    n = `1` + `0` + `<NA>`,
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  )
p1 <- ggplot(tdat, aes(Country, n)) +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Country of enrolment")
p2 <- ggplot(tdat, aes(Country, p_1)) +
  geom_point() +
    labs(
      y = "Proportion\nshortness of breath\nday 28", 
      x = "Country of enrolment")  +
  ylim(0, NA)
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-6-country.pdf"), p, height = 2, width = 6)
p

Site

Relationship site to outcome
tdat <- all_dat %>%
  filter_acs_itt() %>%
  dplyr::count(
    Country_lab = Country,
    Site_lab = fct_infreq(Location),
    Country = factor(PT_CountryName, levels = c("India", "Australia", "Nepal", "New Zealand")),
    Site = PT_LocationName,
    out_sob) %>%
  group_by(Country, Site) %>%
  spread(out_sob, n, fill = 0) %>%
  mutate(
    n = `1` + `0` + `<NA>`,
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  ) %>%
  ungroup()

p1 <- ggplot(tdat, aes(Site_lab, n)) +
  facet_grid( ~ Country, scales = "free_x", space = "free_x") +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        panel.border = element_rect(fill = NA))
p2 <- ggplot(tdat, aes(Site_lab, p_1)) +
  facet_grid( ~ Country, scales = "free_x", space = "free_x") +
  geom_point() +
    labs(
      y = "Proportion\nshortness of breath\nday 28", 
      x = "Site of enrolment") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        panel.border = element_rect(fill = NA))
p <- p1 / p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-6-country-site.pdf"), p, height = 4, width = 6.25)
p

Calendar Time

Relationship calendar date to outcome
caldat <- all_dat %>% 
  filter_acs_itt() %>%
  dplyr::count(out_sob, yr = year(RandDate), mth = month(RandDate)) %>% 
  spread(out_sob, n, fill = 0) %>% 
  mutate(p = `1` / (`1` + `0`),
         n = `1` + `0` + `<NA>`)
p1 <- ggplot(caldat, aes(mth, p)) +
  facet_grid( ~ yr, drop = T, scales = "free_x", space = "free") +
    geom_point() +
    labs(
      y = "Proportion\nshortness of breath\nday 28", 
      x = "Calendar date (month of year)") +
  scale_x_continuous(breaks = 1:12) +
  ylim(0, NA)
p2 <- ggplot(caldat, aes(mth, n)) +
  facet_grid( ~ yr, drop = T, scales = "free_x", space = "free") +
    geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Calendar date (month of year)") +
  scale_x_continuous(breaks = 1:12)
p <- p2 | p1
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-6-calendar-time.pdf"), p, height = 2, width = 6)
p

Modelling

The analysis for shortness of breath is specified equivalently to that for the primary outcome. It includes intervention as randomised, age category, country, site nested within country, epoch, and intervention eligibility. The main analysis restricts the population to only those participants randomised to the anticoagulation domain.

Primary model

Data and sampling
seed <- 36127
mdat <- make_stan_data(
    dat = acs_itt_nona_dat, 
    vars    = c("inelgc3", "agegte60", "ctry"),
    outcome = "out_sob",
    beta_sd = c(10, 2.5, 1, 1))
snk <- capture.output(
  mfit <- model_full$sample(
    data = mdat,
    seed = seed,
    refresh = 0,
    iter_warmup = 1000,
    iter_sampling = 2500,
    chains = 8,
    parallel_chains = 8,
    adapt_delta = 0.95
))
mfit$save_object(file.path("outputs", "models", "secondary", "7-6.rds"))
mdrws <- as_draws_rvars(mfit$draws())
names(mdrws$beta) <- colnames(mdat$X)
names(mdrws$gamma_epoch) <- acs_itt_nona_dat %>% dplyr::count(epoch, epoch_lab) %>% pull(epoch_lab)
names(mdrws$gamma_site) <- acs_itt_nona_dat %>% dplyr::count(site_num, site) %>% pull(site)
site_facet <- factor(mdat$region_by_site, labels = c("India", "Australia\nNew Zealand", "Nepal"))

mdrws$Acon <- attr(mdat$X, "contrasts")$randA %**% mdrws$beta[grepl("randA[0-9]", names(mdrws$beta))]
mdrws$Ccon <- attr(mdat$X, "contrasts")$randC %**% mdrws$beta[grepl("randC[0-9]", names(mdrws$beta))]
mdrws$Atrt <- mdrws$Acon[-1] - mdrws$Acon[1]
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
mdrws$compare <- c(
  "Intermediate vs low" = exp(mdrws$Ctrt[1]),
  "Low with aspirin vs low" = exp(mdrws$Ctrt[2]),
  "Therapeutic vs low" = exp(mdrws$Ctrt[3]),
  "Intermediate vs low with aspirin" = exp(mdrws$Ctrt[1] - mdrws$Ctrt[2]),
  "Intermediate vs therapeutic" = exp(mdrws$Ctrt[1] - mdrws$Ctrt[3]),
  "Low with aspirin vs therapeutic" = exp(mdrws$Ctrt[2] - mdrws$Ctrt[3])
)
mdrws$OR <- c(
  setNames(exp(mdrws$Ctrt), c("Intermediate", "Low with aspirin", "Therapeutic")),
  setNames(exp(mdrws$beta[7:8]), c("Ineligible aspirin", "Age \u2265 60")),
  setNames(exp(mdrws$beta[9:10]), c("Australia/New Zealand", "Nepal"))
)
Odds ratio summary table
save_tex_table(
  odds_ratio_summary_table(mdrws$OR, "latex"),
  "outcomes/secondary/7-6-primary-model-acs-itt-summary-table")
odds_ratio_summary_table(mdrws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Intermediate 0.77 (0.55, 1.09) 0.79 (0.14) 0.93
Low with aspirin 1.19 (0.77, 1.84) 1.22 (0.28) 0.22
Therapeutic 0.93 (0.35, 2.30) 1.02 (0.50) 0.57
Ineligible aspirin 1.19 (0.40, 3.31) 1.36 (0.77) 0.38
Age ≥ 60 2.09 (1.48, 2.94) 2.12 (0.37) 0.00
Australia/New Zealand 2.32 (0.59, 8.49) 2.89 (2.20) 0.11
Nepal 0.52 (0.13, 2.52) 0.72 (0.71) 0.81
Posterior summaries for model parameters (fixed-effects).
Odds ratio summary for epoch and site
p <- plot_epoch_site_terms(
  exp(mdrws$gamma_epoch),
  exp(mdrws$gamma_site),
  site_facet
)
pth <- file.path(
  "outputs", "figures", "outcomes", "secondary", 
  "7-6-acs-itt-epoch-site-terms.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Code
p <- plot_or_densities(mdrws$compare)
p

Code
mfit$summary(variables = c("beta", "gamma_site", "gamma_epoch"))
# A tibble: 48 × 10
   variable    mean  median    sd   mad      q5     q95  rhat ess_bulk ess_tail
   <chr>      <dbl>   <dbl> <dbl> <dbl>   <dbl>   <dbl> <dbl>    <dbl>    <dbl>
 1 beta[1]  -2.01   -2.02   0.611 0.597 -3.00   -0.996   1.00    4875.    8480.
 2 beta[2]  -0.176  -0.170  0.362 0.356 -0.781   0.418   1.00   14046.   14904.
 3 beta[3]   0.0482  0.0470 0.179 0.177 -0.243   0.345   1.00   14419.   15222.
 4 beta[4]  -0.250  -0.251  0.213 0.215 -0.598   0.0984  1.00   13703.   15137.
 5 beta[5]   0.920   0.913  0.539 0.535  0.0405  1.82    1.00   15505.   13871.
 6 beta[6]  -0.129  -0.126  0.320 0.320 -0.659   0.389   1.00   25602.   15513.
 7 beta[7]   0.163   0.172  0.541 0.538 -0.743   1.04    1.00   27357.   15176.
 8 beta[8]   0.736   0.736  0.175 0.176  0.448   1.03    1.00   26772.   14401.
 9 beta[9]   0.833   0.842  0.681 0.673 -0.313   1.93    1.00    7899.   11493.
10 beta[10] -0.631  -0.655  0.749 0.720 -1.82    0.652   1.00   10033.   12915.
# … with 38 more rows
Code
mfit$diagnostic_summary()
$num_divergent
[1] 0 2 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.8071 0.8333 0.7741 0.8229 0.7515 0.7539 0.7825 0.7153
Code
mcmc_trace(mdrws["beta"])

Code
mcmc_trace(mdrws["gamma_site"])

Code
mcmc_trace(mdrws["gamma_epoch"])

Code
mcmc_trace(mdrws[c("tau_site", "tau_epoch")])

Posterior predictive

Code
y_ppc <- mdrws$y_ppc
ppc_dat <- bind_cols(acs_itt_nona_dat, tibble(y_ppc = y_ppc))

grp_ppc <- function(grp) {
  ppc_dat %>%
    group_by(grp = !!grp) %>%
    summarise(
      mean_y = mean(out_sob),
      rvar_mean_y_ppc = rvar_mean(y_ppc)
    )
}
plot_grp_ppc <- function(dat, lab = "") {
  dat %>%
    ggplot(aes(y = grp, xdist = rvar_mean_y_ppc)) +
    stat_interval(size = 2) +
    geom_point(aes(x = mean_y, y = 1:nrow(dat)), colour = "red", shape = 23) +
    labs(
      x = "Posterior predictive\nproportion", 
      y = lab)  
}

ppc_C <- grp_ppc(sym("CAssignment"))
ppc_ctry <- grp_ppc(sym("country"))
ppc_epoch <- grp_ppc(sym("epoch"))
ppc_site <- ppc_dat %>%
  group_by(Country = country, grp = site_raw) %>%
  summarise(
    mean_y = mean(out_sob),
    rvar_mean_y_ppc = rvar_mean(y_ppc)
  )

p1 <- plot_grp_ppc(ppc_C, "Anticoagulation\nintervention") + labs(x = "")
p2 <- plot_grp_ppc(ppc_ctry, "Country") + labs(x = "")
p3 <- plot_grp_ppc(ppc_epoch, "Epoch") + labs(x = "")
p4 <- plot_grp_ppc(ppc_site %>% filter(Country == "IN"), "Sites\nIndia")  + labs(x = "")
p5 <- plot_grp_ppc(ppc_site %>% filter(Country == "AU"), "Sites\nAustralia") + labs(x = "")
p6 <- plot_grp_ppc(ppc_site %>% filter(Country == "NP"), "Sites\nNepal")
p7 <- plot_grp_ppc(ppc_site %>% filter(Country == "NZ"), "Sites\nNZ")
p <- ((p3 | p1 / p2) / 
        ( (p4 | p5) / (p6 | p7) + 
            plot_layout(heights = c(3, 1)) ) ) +
  plot_layout(heights = c(1, 1.5), guides = "collect")
pth <- file.path("outputs", "figures", "outcomes", "secondary",
                 "7-6-primary-model-acs-itt-ppc.pdf")
ggsave(pth, p, width = 6, height = 5.5)
p

Sensitivity: Concurrent Controls

A sensitivity analysis is performed on the reduced concurrent control analysis sets for each intervention.

Intermediate

Overview of outcome
save_tex_table(
  make_summary_table(acs_itt_concurc2_dat, "latex"),
  "outcomes/secondary/7-6-anticoagulation-concurrent-intermediate-summary")
make_summary_table(acs_itt_concurc2_dat)
n (\%) Patients Known Missing Shortness of breath day 28
Low dose 610 577 (94.6) 33 (5.4) 115 (19.9)
Intermediate dose 613 584 (95.3) 29 (4.7) 110 (18.8)
Overall 1223 1161 (94.9) 62 (5.1) 225 (19.4)
Summary of outcome for intermediate dose concurrent enrolments.
Fit model
mdat <- make_stan_data_concurrent(
  acs_itt_concurc2_nona_dat,
  c("agegte60", "ctry"),
  "out_sob",
  c(2.5, 1, 1)
)
snk <- capture.output(
  mfit <- model_fixed_only$sample(
    data = mdat,
    seed = 38135,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500,
    refresh = 0
))
mdrws <- as_draws_rvars(mfit$draws(c("beta")))
names(mdrws$beta) <- colnames(mdat$X)
mdrws$Ccon <- attr(mdat$X, "contrasts")$randC %**% mdrws$beta[2]
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
mdrws$compare <- c("Intermediate vs low" = exp(mdrws$Ctrt[1]))
mdrws$OR <- c(
  setNames(exp(mdrws$Ctrt), c("Intermediate")),
  setNames(exp(mdrws$beta[-(1:2)]), c("Age \u2265 60", "Australia/New Zealand", "Nepal"))
)
Odds ratio summary
save_tex_table(
  odds_ratio_summary_table(mdrws$OR, "latex"),
  "outcomes/secondary/7-6-primary-model-acs-itt-concurrent-intermediate-summary-table")
odds_ratio_summary_table(mdrws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Intermediate 0.91 (0.67, 1.23) 0.92 (0.14) 0.74
Age ≥ 60 1.64 (1.19, 2.24) 1.66 (0.27) 0.00
Australia/New Zealand 3.91 (2.62, 5.82) 4.00 (0.82) 0.00
Nepal 0.47 (0.22, 0.92) 0.50 (0.18) 0.99
Odds ratio summary table for model.
Posterior contrast
p <- plot_or_densities(mdrws$compare)
p

Code
mfit$summary()
# A tibble: 11 × 10
   variable         mean   median     sd    mad       q5      q95  rhat ess_bulk
   <chr>           <dbl>    <dbl>  <dbl>  <dbl>    <dbl>    <dbl> <dbl>    <dbl>
 1 lp__        -550.     -5.49e+2 1.57   1.38   -5.53e+2 -548.     1.00    9385.
 2 beta_raw[1]   -0.689  -6.88e-1 0.0401 0.0403 -7.55e-1   -0.623  1.00   17830.
 3 beta_raw[2]   -0.0687 -6.95e-2 0.107  0.107  -2.44e-1    0.107  1.00   23339.
 4 beta_raw[3]    0.198   1.99e-1 0.0645 0.0648  9.05e-2    0.303  1.00   19354.
 5 beta_raw[4]    1.37    1.36e+0 0.203  0.199   1.03e+0    1.70   1.00   21689.
 6 beta_raw[5]   -0.758  -7.46e-1 0.364  0.359  -1.38e+0   -0.182  1.00   25118.
 7 beta[1]       -1.72   -1.72e+0 0.100  0.101  -1.89e+0   -1.56   1.00   17830.
 8 beta[2]       -0.0687 -6.95e-2 0.107  0.107  -2.44e-1    0.107  1.00   23339.
 9 beta[3]        0.495   4.96e-1 0.161  0.162   2.26e-1    0.756  1.00   19354.
10 beta[4]        1.37    1.36e+0 0.203  0.199   1.03e+0    1.70   1.00   21689.
11 beta[5]       -0.758  -7.46e-1 0.364  0.359  -1.38e+0   -0.182  1.00   25118.
# … with 1 more variable: ess_tail <dbl>
Code
mfit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 1.049 1.023 1.036 1.120 1.035 1.026 1.071 1.012
Code
mcmc_trace(mdrws["beta"])

Low-dose with aspirin

Overview of outcome
save_tex_table(
  make_summary_table(acs_itt_concurc3_dat, "latex"),
  "outcomes/secondary/7-6-anticoagulation-concurrent-stdaspirin-summary")
make_summary_table(acs_itt_concurc3_dat)
n (\%) Patients Known Missing Shortness of breath day 28
Low dose 299 280 (93.6) 19 (6.4) 55 (19.6)
Intermediate dose 298 280 (94.0) 18 (6.0) 52 (18.6)
Low dose with aspirin 283 271 (95.8) 12 (4.2) 59 (21.8)
Overall 880 831 (94.4) 49 (5.6) 166 (20.0)
Summary of outcome for low-dose with aspirin concurrent enrolments.
Fit model
mdat <- make_stan_data_concurrent(
  acs_itt_concurc3_nona_dat,
  c("agegte60", "ctry"),
  "out_sob",
  c(2.5, 1)
)
snk <- capture.output(
  mfit <- model_fixed_only$sample(
    data = mdat,
    seed = 13578,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500,
    refresh = 0
))
mdrws <- as_draws_rvars(mfit$draws(c("beta")))
names(mdrws$beta) <- colnames(mdat$X)
mdrws$Ccon <- attr(mdat$X, "contrasts")$randC %**% mdrws$beta[2:3]
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
mdrws$compare <- c("Intermediate vs low" = exp(mdrws$Ctrt[1]),
                   "Low with aspirin vs low" = exp(mdrws$Ctrt[2]),
                   "Intermediate vs low with aspirin" = exp(mdrws$Ctrt[1] - mdrws$Ctrt[2]))
mdrws$OR <- c(
  setNames(exp(mdrws$Ctrt), c("Intermediate", "Low with aspirin")),
  setNames(exp(mdrws$beta[-(1:3)]), c("Age \u2265 60", "Australia/New Zealand"))
)
Odds ratio summary
save_tex_table(
  odds_ratio_summary_table(mdrws$OR, "latex"),
  "outcomes/secondary/7-6-primary-model-acs-itt-concurrent-stdaspirin-summary-table")
odds_ratio_summary_table(mdrws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Intermediate 0.96 (0.64, 1.46) 0.99 (0.21) 0.57
Low with aspirin 1.17 (0.78, 1.77) 1.20 (0.26) 0.22
Age ≥ 60 1.54 (1.05, 2.23) 1.56 (0.30) 0.01
Australia/New Zealand 2.17 (1.17, 3.91) 2.26 (0.70) 0.01
Odds ratio summary table for model.
Posterior contrast
p <- plot_or_densities(mdrws$compare)
p

Code
mfit$summary()
# A tibble: 11 × 10
   variable         mean   median     sd    mad       q5      q95  rhat ess_bulk
   <chr>           <dbl>    <dbl>  <dbl>  <dbl>    <dbl>    <dbl> <dbl>    <dbl>
 1 lp__        -417.     -4.17e+2 1.58   1.42   -4.20e+2 -415.     1.00    9210.
 2 beta_raw[1]   -0.631  -6.30e-1 0.0447 0.0450 -7.05e-1   -0.558  1.00   16753.
 3 beta_raw[2]    0.140   1.42e-1 0.150  0.151  -1.08e-1    0.386  1.00   22687.
 4 beta_raw[3]   -0.0510 -4.97e-2 0.150  0.152  -2.99e-1    0.195  1.00   21260.
 5 beta_raw[4]    0.170   1.71e-1 0.0770 0.0764  4.16e-2    0.296  1.00   18995.
 6 beta_raw[5]    0.768   7.74e-1 0.307  0.303   2.55e-1    1.27   1.00   19666.
 7 beta[1]       -1.58   -1.57e+0 0.112  0.112  -1.76e+0   -1.39   1.00   16753.
 8 beta[2]        0.140   1.42e-1 0.150  0.151  -1.08e-1    0.386  1.00   22687.
 9 beta[3]       -0.0510 -4.97e-2 0.150  0.152  -2.99e-1    0.195  1.00   21260.
10 beta[4]        0.426   4.29e-1 0.193  0.191   1.04e-1    0.740  1.00   18995.
11 beta[5]        0.768   7.74e-1 0.307  0.303   2.55e-1    1.27   1.00   19666.
# … with 1 more variable: ess_tail <dbl>
Code
mfit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 1.127 1.112 1.099 1.106 1.162 1.131 1.058 1.078
Code
mcmc_trace(mdrws["beta"])

Therapeutic

Overview of outcome
save_tex_table(
  make_summary_table(acs_itt_concurc4_dat, "latex"),
  "outcomes/secondary/7-6-anticoagulation-concurrent-therapeutic-summary")
make_summary_table(acs_itt_concurc4_dat)
n (\%) Patients Known Missing Shortness of breath day 28
Low dose 79 72 (91.1) 7 (8.9) 18 (25.0)
Intermediate dose 65 61 (93.8) 4 (6.2) 13 (21.3)
Therapeutic dose 50 44 (88.0) 6 (12.0) 11 (25.0)
Overall 194 177 (91.2) 17 (8.8) 42 (23.7)
Summary of outcome for therapeutic dose concurrent enrolments.
Fit model
mdat <- make_stan_data_concurrent(
  acs_itt_concurc4_nona_dat,
  c("agegte60", "ctry"),
  "out_sob",
  c(2.5, 1, 1)
)
snk <- capture.output(
  mfit <- model_fixed_only$sample(
    data = mdat,
    seed = 13578,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500,
    refresh = 0
))
mdrws <- as_draws_rvars(mfit$draws(c("beta")))
names(mdrws$beta) <- colnames(mdat$X)
mdrws$Ccon <- attr(mdat$X, "contrasts")$randC %**% mdrws$beta[2:3]
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
mdrws$compare <- c("Intermediate vs low" = exp(mdrws$Ctrt[1]),
                   "Therapeutic vs low" = exp(mdrws$Ctrt[2]),
                   "Intermediate vs therapeutic" = exp(mdrws$Ctrt[1] - mdrws$Ctrt[2]))
mdrws$OR <- c(
  setNames(exp(mdrws$Ctrt), c("Intermediate", "Therapeutic")),
  setNames(exp(mdrws$beta[-(1:3)]), c("Age \u2265 60", "India", "Australia/New Zealand"))
)
Odds ratio summary
save_tex_table(
  odds_ratio_summary_table(mdrws$OR, "latex"),
  "outcomes/secondary/7-6-primary-model-acs-itt-concurrent-therapeutic-summary-table")
odds_ratio_summary_table(mdrws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Intermediate 0.46 (0.18, 1.13) 0.51 (0.25) 0.95
Therapeutic 0.66 (0.25, 1.68) 0.74 (0.38) 0.80
Age ≥ 60 1.89 (0.84, 4.34) 2.07 (0.93) 0.06
India 1.25 (0.36, 4.03) 1.49 (0.96) 0.36
Australia/New Zealand 12.28 (5.66, 28.30) 13.51 (5.93) 0.00
Odds ratio summary table for model.
Posterior contrast
p <- plot_or_densities(mdrws$compare)
p

Code
mfit$summary()
# A tibble: 13 × 10
   variable    mean  median    sd   mad       q5     q95  rhat ess_bulk ess_tail
   <chr>      <dbl>   <dbl> <dbl> <dbl>    <dbl>   <dbl> <dbl>    <dbl>    <dbl>
 1 lp__     -83.2   -82.9   1.76  1.60  -86.6    -81.0    1.00    9745.   12606.
 2 beta_ra…  -1.03   -1.03  0.149 0.148  -1.29    -0.797  1.00   10110.   11479.
 3 beta_ra…   0.256   0.259 0.355 0.358  -0.333    0.838  1.00   19380.   14705.
 4 beta_ra…   0.490   0.488 0.332 0.327  -0.0542   1.05   1.00   17811.   14955.
 5 beta_ra…   0.256   0.255 0.168 0.169  -0.0197   0.533  1.00   13647.   12806.
 6 beta_ra…   0.214   0.225 0.615 0.616  -0.809    1.21   1.00   16860.   14995.
 7 beta_ra…   2.52    2.51  0.409 0.408   1.86     3.20   1.00   12324.   12874.
 8 beta[1]   -2.58   -2.57  0.373 0.371  -3.22    -1.99   1.00   10110.   11479.
 9 beta[2]    0.256   0.259 0.355 0.358  -0.333    0.838  1.00   19380.   14705.
10 beta[3]    0.490   0.488 0.332 0.327  -0.0542   1.05   1.00   17811.   14955.
11 beta[4]    0.640   0.638 0.420 0.422  -0.0493   1.33   1.00   13647.   12806.
12 beta[5]    0.214   0.225 0.615 0.616  -0.809    1.21   1.00   16860.   14995.
13 beta[6]    2.52    2.51  0.409 0.408   1.86     3.20   1.00   12324.   12874.
Code
mfit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 1.026 1.017 1.050 1.108 1.078 1.017 1.064 1.034
Code
mcmc_trace(mdrws["beta"])